What is the Entanglement Length in a Polymer Melt? 
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We present the results of molecular dynamics simulations of very long model polymer chains analyzed by 
various experimentally relevant techniques. The segment motion of the chains is found to be in very good 
agreement with the reptation model. We also calculated the plateau modulus G%. The predictions of the 
entanglement length A^e from G% and from the mean square dispacement of the chain segments disagree by a 
factor of about 2.2(2), indicating an error in the prefactor in the standard formula for G%. We show that recent 
neutron spin echo measurements were carried out for chain lengths which are too small to allow for a correct 
determination of A^e. 



How an entangled polymer chain moves in a dense melt 
of other chains has been a long standing problem. The most 
widely accepted picture is that the chains reptate like a snake 
[|l|,|[l. For short times chains move isotropically until they feel 
the constraints of their neighboring chains. For intermediate 
times, the chain segments move along the path or tube created 
by the surrounding chains in a Rouse-like ^] motion. Only 
the ends explore new space. For the inner section of the chain, 
this Rouse-like motion on a contour which is also a random 
walk gives rise to the famous t^^'^ power law regime in the 
mean-square displacement of the beads gi{t). For very long 
time, the motion is diffusive, with a chain diffusion constant 
D that scales as N^'^ for large N, where N is the chain length. 
For short chains, the motion is much simpler and can be ap- 
proximately described by the Rouse model and D ~ N^^. 
To characterize the crossover between the Rouse and reptation 
regime, one can define an entanglement length Ne- Within the 
reptation model, N^ can be related to both the tube diameter 
dr and crossover time t^ from the early time Rouse regime 
where gi{t) ~ t^^^ to the t^/"* regime as well as from the 
value of the plateau modulus G^. Over the years there have 
been a number of experiments [H-h] and simulations [n-Jntl 
designed to test various aspects of the theory. Recent neutron 
spin echo (NSE) scattering experiments [^] which measure 
the dynamic structure factor S{k, t), suggest that N^ as mea- 
sured from (It is consistent with that determined from G'j^. 
Our previous simulation results [0,^ suggested an inconsis- 
tency in that Ne measured from gi (i) was about one half that 
measured from S{k, t), though both gi{t) and S{k, t) are sin- 
gle chain quantities and measure the same motion. However 
since the chains were only a few N^ (N < 200), our results 
were not conclusive. We have now extended our simulations 
to much longer chains, as long as A^ = 10000, and measured 
not only single chain quantities but also G^. We find clear ev- 
idence for differences on the order two between N^ measured 
from gi{t) or S(k,t) and that from G^. The previous re- 
ported agreement in N^ determined from NSE data for S{k, t) 
and G% is shown to be largely due to finite chain length ef- 
fects in dr as determined from S{k,t). 

To overcome the long time scales needed to simulate a melt 
of long entangled polymers, we use a coarse grained model 



in which the polymer is treated as a string of beads of mass 
m connected by a spring. The beads interact with pure re- 
pulsive Lennard-Jones excluded volume interactions (cutoff 
at 2^/^ct) and are connected by a finite extensible non-linear 
elastic potential (FENE) between neighbors along the chain 
(see e.g. [Pfi for details). The model parameters are the same 
as in ref. [0]. The temperature T = e/ks, where e is the 
strength of the Lennard-Jones interaction. We use dimension- 
less units in which a = 1 and e = 1 and the basic unit of time 

We performed constant volume simulations of monodis- 
perse polymer melts at a segment density of p = 0.85cr^^. 
The temperature was kept constant by coupling the motion of 
the beads weakly to a heat bath with a local friction coefficient 
r = 0.5t^^. The equations of motion were integrated using 
a velocity Verlet algorithm with a time step At = 0.012t. 
The average bond distance is \/(p) ~ 0.97a and chain stiff- 
ness Coo — 1.75. This gives a statistical segment length of 
b = 1.28ct. Initial conformations of the chains were grown 
as non-reversal random walks with the proper melt end-to- 
end extension. Resulting inital overlaps of chain segments 
were removed by simulating a soft core potential for a very 
short time to avoid instabilities with the Lennard-Jones poten- 
tial. Initially [0], we studied chains of length 5 < iV < 200 
and later reran [^ all the systems for N < 200 with more 
chains for longer times to improve the quality of the data. Our 
new results are for a system of M chains of length N, for 
M/N = 120/350, 350/700 and 50/10000. 

The most direct route to verify the predictions of the rep- 
tation model is to monitor the mean-square displacements of 
the segments ri, 

ffz,i(i) = ((r.(i)-r.(0))') (1) 

g^At) = ((r,(t) - r„„(t) - r,(0) + r,„(0))') , (2) 



and the center of mass of the chains Vcm (t), 
gsit) = ({remit) - rcm{0)f) 



(3) 



The reptation model predicts the following power laws for var- 
ious time regimes [| 
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FIG. 1. Mean-square displacements 32 (t) (open symbols) and 
g3{t) (closed symbols) for chain length TV = 350 (□), N = 700 (o) 
and N = 10000 (A). The straight lines show some power law be- 
haviors to guide the eye. The local reptation power laws 32 (i) oc t^''^ 
and 33(4) oc t^'^ are verified with remarkable clarity. 
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where W = ^^n^^ dr the effective tube diameter and C is 
the effective bead friction. 52 (i) shows the same regimes for 
t < Td, but goes to a plateau value of Rq{N) for t > Td- 
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where Te is the entanglement time, tr the Rouse time and Td 
the disentanglement time [0]: 
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In Fig. |l] we show the results for 52 (i) for the innermost seg- 
ments of the chains and (73 (t) for the systems with chain length 
N = 350, 700 and 10000. While the simulation times for 
N — 350 were long enough to reach the diffusive regime the 
data for N — 700 and 10000 just reach far into the predicted 
reptation regime. 

After an initial Rouse-like motion 52 (i) oc t^/'^ for chains 
segments up to a time Te = 1420 ± IOOt {N = 1100 ± IOOt) 
for N = 700 {N = 10000) the motion slows down and is pro- 
portional to t°^^'^' which is in remarkable agreement to the 
reptation model and is in less accordance with Schweizer's 
mode coupling theory [n2p . The crossover time Te leads us to 
our first estimate for A^e by assuming it to be the Rouse relax- 
ation time of a subchain of length A^e- From the initial slope 



of the short time Rouse-regime, gi{t) = 0.525(5)cr^(i/T)^/^, 
we determine W = 0.025 (2)Tcr~^ which is consistent with a 
bead friction of C = 25(1)t^^ determined separately by the 
relaxation of the Rouse modes of shorter chains | ]13[ | . Insert- 
ing this is into the expression for tr, Eq. (@), yields A^e = 
32(2) for N = 700 and A^e = 28(2) for N = 10000. By 
equating the initial two power-law regimes for gi{t) at Te, we 
determine a tube diameter dx = 7.6(3) for N — 700. Note, 
that for A^ ~ 10000 the prefactor of the i°'^^ regime appears 
about 7% smaller which gives dr = 7.1(3). Assuming that 
dr = R'^{Ne) [0], where R^ is the end-to-end distance, gives 
A^e = 35(2) for N = 700 and N^ = 32(2) for N = 10000. 
Thus the two ways of defining A"e give consistent results, d^ 
is also proporptional to gi(Te), though the exact prefactor is 
not strictly specified. Employing a Gaussian picture [^ of the 
tube one can estimate 51 (te) = 2Rc{N^) — d\,/i. With the 
values of 51 (Te, 700) = 18.9(5) and 51 (Te, 10000) = 17.4(5) 
one obtains A^e = 35(1), dr = 7.5(2) and A^e = 32(l),dT = 
7.2(2) repsectively. These values agree with our old results 
[Q] within eiTor bars. After about the Rouse time tr{N) the 
dynamics of g2{t) should cross over to a second t^ regime, 
which corresponds to the diffusion of the whole chain along 
the gaussian tube contour. This second regime is not visible 
for A^ = 350 since the chains are not long enough and only a 
broad crossover to the final plateau is observed. This regime 
should be more pronounced for A^ = 700, but the computa- 
tional effort to obtain it is prohibitively large at present (about 
a CPU month on a 256 processor T3E). The slightly subdif- 
fusive behavior of 53 (i) for times shorter than Te is not due 
to entanglement effects and will be discussed elsewhere Jl3|]. 
After Te a clear f^/^ regime in gsit) is observed for A^ — 700 
and 10000, in agreement with the reptation model rather than 
mode coupling [O]. The ratio of the power-law prefactors for 
these two chainlengths is 18(2). This is in good agreement 
again with the reptation model Eq. (0), where we expect a 
ratio 16.3 taking the slight A^-dependence of dx into account. 
32 (i) for the shorter A^ — 350 chains show a slightly higher 
exponent of t°-^'^i'^)_ After about 3.5 • IO^t, about twice the 
Rouse time Tfl(350) — 1.8 • IO'^t, the data show diffusive 
behavior. 

Experimentally the motion of the segments can be obtained 
by measuring the time-dependent single-chain structure func- 
tion 



^(J''^) = j^ /^exp(*k. (r,(t) -r,(0)))\ . (7) 

\ i,j I 

For reptating chains this is predicted to be of the approxima- 
;ive forn 
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tive form in the limits -|^ < fc < 4^ and t > t^: 
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where f{u) — exp(u^/36)erfc(M/6). The short time. Rouse- 
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FIG. 2. Dynamic single chain structure function S{k,t) for dif- 
ferent chain length TV = 350 (+), 700 (x), and and 10000 (■) for 
various fc-values. The solid lines are fits to Eq. (tel) (dr given in the 
tex). For equal fc-values the plateaus show a strong A'^ dependence. 
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FIG. 3. Normal tension ajv measured as a function of time after 
a step-train for various strain amplitudes A (1.25 (□), 1.5 (o), 1.75 
(A) and 2.0 (O)). The straight lines are fitted by hand to the long 
time decay of the stress to extrapolate the plateau stress. 



like motion is not described by this formula and only en- 
ters through the inverse friction coefficient W. The origi- 
nal formula was derived by de Gennes [|l^ on the basis of 
the reptation model and was slightly modified by Schlegel 
et al. [^ in an attempt to extend the range of validity of 
the original expression to larger fc-vectors. Eq. (g) dif- 
fers somewhat from the form given in ref. [||, as its time- 
dependence has been corrected due an argument by Kremer 
and Binder [|l5|]. For the present set of data this correction 
is very small and can be neglected. Since for our system, 
dr — 7a, the upper bound in fc-space kmax — 1-Ocr^^. The 
lower bound for each N with Rq — b'^N/6 is kmin{N) ~ 
0.1(10000), 0.4(700), 0.8(350)cr-i. 

If one inserts r^, Eq. (^, in Eq. (Q) the resulting expression 
for S{k, t) contains only one adjustable parameter, dr- We 
calculated S{k, t) for N = 350, 700 and 10000 for several k- 
values between 0.2 < fccr < 1.0 and fitted the data to Eq. (ph 
in the time window 5000r < i < lOOOOOr. The best fit gives 
dr = 15.7(5)cr for N = 350 (fc = QAa-^), dr = 12.7(3)<t 
for N = 700 (fc == 0.4cr-i) ^^^ ^^ ^ 8.5(3)cr for N = 
10000 in a simultaneous fit to fc = 0.2, 0.3and0.4cr^^. Figure 
shows our results together with the fitting curves. One can 
see that the agreement of Eq. (pi) is only acceptable for N — 
10000 and N — 700 and is rather poor for the shorter chains. 
The large difference between the tube diameter obtained from 
the data for different chain length suggests that finite chain 
length effects are much more important for S{k,t) than in 
gi{t) in determining dr- These finite chain length effects are 
not accounted for by Eq. (g). Clearly, the apparent value of 
dx approaches our previous estimates of dx, which should be 
expected since both methods measure the same quantity. For 
N = 10000 finite chain length effects should be very small. 
Assuming a finite size scaling of d{N) — doo+a/Ny a simple 
fit gives doc = 7.65 and y = 0.67 showing that finite size 



effects decay very slowly and the extrapolated estimate for dx 
aggrees nicely with our estimate from gi (re). 

The standard method to determine 7Ve,p (for clarity we in- 
dex Ne determined from the plateau-modules with an addi- 
tional index p) experimentally is by measuring the plateau 
modulus in an oscillatory shear experiment. Alternatively, it 
is also possible to measure the normal stress decay in a step 
strain elongation. Since the latter is much simpler to perform 
in a simulation we ran volume conserving step strains for four 
different amplitudes A = 1.25, 1.5, 1.75 and 2.0. After a rapid 
decay at short times, the stress had a well defined plateau from 
which we could determine G% (see Fig. ^. The normal stress 
o'N — cra;a; — ((7j,y+cr22)/2 was determined by the microscopic 
virial-tensor, x being the direction of elongation. 

We fitted our results to the stress-strain formulas for clas- 
sical rubber elasticity Jig ] ctn = G (A^ — j) and to the 
Mooney-Rivlin (MR) formula ctat = 2Gi (A^ - j) + 
2G2 (A - J3-) to determine G%. The fit to the the MR for- 
mula is excellent and gives G^ = 0.0105kBTa~^ while the 
classical fit is fair and gives a value of G% — O.OOSfcsTa"^. 
It is known experimentally that the MR formula slightly over- 
estimates the modulus while the classical equation always un- 
derestimates it. The standard formula of Doi [p|Jl8|] to calcu- 
late N, 
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gives Ne^p = 65 for the MR fit and Ne,p = 80 for the clas- 
sical formula. Both values are much higher than our previous 
estimate, N^ = 32. 

If one scales the diffusion constant D{N) by the Rouse 
diffusion constant Dji{N) and plots it versus N/N^^p ex- 
perimental results for different polymers [p^[-pl|| and simu- 
lation results for different models [p[pl]] fall onto the same 
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FIG. 4. Scaled diffusion constant D{N)/Dr{N) vs. scaled 
cliain lengtli N/Ne,p for polystyrene (•) ^^ {Me,p = 14600, 
T = 485is:), polyethylene (■) ||9| (M^,p = 870, T = 448ii:), 
PEB2 (a) [0 (M^,p = 992, T = 448^"), our bead spring model 
(A) (Ne,p = 72), the bond-fluctuation model for ^ = 0.5 (□) [|] 
and tangent hard spheres at $ = 0.45 (o) [pl|. All data are scaled 
with A^e^p from the plateau modulus or with 2.27Ve from gi (t). 



ment to the predictons of the reptation model. The dynamical 
exponent of t^/'' for the local reptation regime has been veri- 
fied with remarkable clarity. We further demonstrate that very 
long chains with N > lOOiVe.p are needed for S{k, t) to ar- 
rive at a consistent prediction of A^e with that from the mean- 
square displacements. The most recent experiments were per- 
formed for chain lengths well below this threshold. In contrast 
the formula for the modulus by Doi leads to an estimate of N^ 
larger by a factor of about 2.2(2). Whether this discrepancy 
is due to just uncertainties in prefactors of the reptation model 
or due to the failure of the classical single chain picture for the 
viscoelasticity still remains unclear. 

Most of the simulations were carried out at the Rechenzen- 
trum of the MPG in Munich and at Exxon Research and En- 
gineering Company. Sandia is a multiprogram laboratory op- 
erated by Sandia Corporation, a Lockheed Martin Company, 
for the United States Department of Energy under Contract 
DE-AC04-94AL85000. 



universal curve (see Fig. Mj WM). The diffusion constants 
for the simulated systems are determined by extrapolating 
(73 (i — > oo). Using a value N^^p — 72, intermediate be- 
tween the MR and classical fits to the stress-strain data, to 
normalize N , our results nicely fall onto the experimental val- 
ues. For the bond-fluctuation [g] and tangent hard sphere 
models Jll] ] no plateau-moduli are available, thus we scaled 
Ne calculated from (?i(Te) in these models (bond-fluctuation: 
Ne = 30, tangent hard sphere: N^ — 29) by the same factor 
of Ngp/N^ = 72/32 = 2.25 obtained from our model to esti- 
mate iVg p. However, due to the uncertainty and limited range 
of the simulation data and the scatter in the experimental data, 
Ne^p in the range 2.0 — 2.4:Ne could also be chosen to collapse 
the data. The nice collapse of the data supports our assump- 
tion that Ne and Ne.p are related by a universal multiplicative 
factor of about 2.2(2). It remains unclear though whether this 
prefactor of about 2.2 is truly universal or just a consequence 
of the fact that all three model systems are fully flexible and 
have almost the same packing fractions. 

In the light of our simulation results one should critically 
review the results of recent NSE experiments [ra] which claim 
to support the reptation prediction and rule out other theories 
by fitting the data to Eq. dq). They also claim that their esti- 
mated value of Ne agrees nicely with the value derived from 
the plateau modulus. However, it should be noted that the 
chain lengths in these investigations are only about twice our 
N = 700 chains, i.e. N w 23iVe,p, in a comparable range of 
/c-vectors. The simulation results suggest, that dr determined 
in these experiments is systematically too high by about a fac- 
tor of 1.5, giving a factor of 2 for N^- Note, that the finite 
chain length effects are much stronger in S{k,t) than in 52 (i)- 

To conclude, we find that our data are in very good agree- 
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